library(tidyverse)
library(sf)
library(kableExtra)
library(units)
library(gt)
library(patchwork)
library(purrr)
library(ggtext)

set.seed(123)

distrito <- read_sf("dados/distrito/SIRGAS_SHP_distrito.shp") %>% st_set_crs("epsg:31983")

Microdados do Censo de 2022

É importante destacar que os resultados disponíveis são parciais do censo, com dados atualizados até o momento da divulgação. No entanto, mesmo parciais, esses dados podem ser extremamente úteis para suas análises.

Um ponto crucial a considerar é que o nível mínimo de observação georreferenciada nos microdados do censo são os setores censitários. Os setores censitários são unidades geográficas definidas pelo Instituto Brasileiro de Geografia e Estatística (IBGE) para coleta e tabulação de dados censitários. Eles são delimitados de forma a garantir uma cobertura completa e homogênea de todo o território nacional, facilitando a análise comparativa entre diferentes áreas geográficas. A delimitação geográfica do setor também considera questões logísticas para o entrevistador conseguir entrevistar a todos.

# Read dados do censo 2022
censo <- read_sf("dados/censo/SP_Malha_Preliminar_2022.shp") %>% 
  filter(CD_MUN == "3550308") %>% 
  st_transform("epsg:31983") %>% # Sistema de coordenadas do geosampa
  select(id_setor_censitario = CD_SETOR, v0001:v0007) %>% 
  mutate(area_setor = st_area(geometry))
id_setor_censitario v0001 v0002 v0003 v0004 v0005 v0006 v0007 geometry area_setor
355030836000023P 523 200 200 0 2.827027 49.189189 185 MULTIPOLYGON (((358863.3 74… 26327.781 [m^2]
355030825000026P 547 200 200 0 2.878947 18.421053 190 MULTIPOLYGON (((356344.7 73… 7180.695 [m^2]
355030833000279P 799 252 252 0 3.400000 11.914894 235 MULTIPOLYGON (((353893.9 73… 46998.556 [m^2]
355030843000013P 499 206 206 0 2.851429 21.714286 175 MULTIPOLYGON (((318788.2 73… 36182.414 [m^2]
355030830000764P 654 238 238 0 3.099526 9.952607 211 MULTIPOLYGON (((328927.3 73… 33462.349 [m^2]
355030890000334P 280 138 138 0 2.258065 20.161290 124 MULTIPOLYGON (((333870.3 73… 13169.581 [m^2]
355030838000607P 188 102 102 0 1.938144 9.278351 97 MULTIPOLYGON (((332865.9 73… 1972.884 [m^2]
355030853000064P 511 302 302 0 2.147059 1.260504 238 MULTIPOLYGON (((338008.3 73… 44051.010 [m^2]
355030809000067P 264 140 138 2 2.357143 17.857143 112 MULTIPOLYGON (((333077.1 73… 26519.305 [m^2]
355030836000443P 793 300 300 0 3.223577 4.065041 246 MULTIPOLYGON (((358244.3 74… 10631.632 [m^2]
gg <- ggplot() +
  geom_sf(data = censo %>%
            mutate(densidade = v0001/area_setor,
                   decil = ntile(densidade, 10)),
          aes(fill = factor(decil)), color = NA) +
  scale_fill_viridis_d() +
  labs(fill = "Decil de densidade") +
  theme_void()

ggsave("tex/imagens/mapa.png", gg, width = 20, height = 20, dpi = 400)

descritiva <- censo %>% 
  st_drop_geometry() %>% 
  summarize("Total de pessoas" = sum(v0001),
            "Total de Domicílios" = sum(v0002),
            "Total de Domicílios Particulares Ocupados" = sum(v0007)) %>% 
  pivot_longer(everything())
name value
Total de pessoas 11.451.999
Total de Domicílios 4.996.529
Total de Domicílios Particulares Ocupados 4.316.336
censo %>% 
  mutate(distancia_centro = st_distance(geometry, 
                                        read_sf("dados/distrito/SIRGAS_SHP_distrito.shp") %>% 
                                          st_set_crs("epsg:31983") %>% 
                                          filter(ds_nome == "SE")) %>% 
           as.numeric() %>% 
           cut(breaks = 10^3*(0:40), labels = FALSE)) %>% 
  st_drop_geometry() %>%
  group_by(distancia_centro) %>% 
  summarize(densidade = sum(v0001) / sum(as.numeric(area_setor) / 10^6)) %>% 
  ggplot(aes(x = distancia_centro, y = densidade)) +
  geom_col() +
  scale_y_continuous(labels = scales::comma_format(big.mark = ".")) +
  scale_x_continuous(labels = scales::comma_format(big.mark = ".")) +
  labs(x = "Distância ao centro em quilômetros", 
       y = "Densidade populacional por quilômetro quadrado") +
  theme_classic()

Base de Dados do IPTU

A base de dados do IPTU (Imposto Predial e Territorial Urbano) é uma fonte abrangente de informações sobre imóveis urbanos dentro do município. Essa base é considerada completa, pois abrange todos os imóveis sujeitos à tributação do IPTU, representando uma fonte confiável e abrangente de informações sobre a propriedade urbana. Propriedades que não foram construídas dentro da legalidade não constam nessa base.

O número do contribuinte, utilizado para identificar exclusivamente cada imóvel na base de dados do IPTU, é representado diretamente pelo SQL, sendo essencial para consultas e manipulação dos dados relacionados ao imposto predial e territorial urbano. Esse formato permite a integração e cruzamento com outras bases de dados, como a base de lotes, que está georreferenciada, fornecendo informações espaciais adicionais que não estão disponíveis na base do IPTU.

Procedimento para Lotes de Condomínios

Quando um lote é um condomínio, ou seja, quando não se classifica de acordo com condominio == "00-0", é necessário substituir os múltiplos números de SQL pelo número do condomínio. Isso ocorre porque cada unidade dentro do condomínio pode ser tratada como uma entidade separada para fins tributários, mas para a base de lotes, estão todos juntos.

IPTU <- read.csv("dados/IPTU/IPTU_2024.csv", sep=";", encoding = "latin1") %>% 
    as_tibble() %>% 
    select(sql = "NUMERO.DO.CONTRIBUINTE", 
           condominio = "NUMERO.DO.CONDOMINIO",
           area_terreno = "AREA.DO.TERRENO",
           area_construida = "AREA.CONSTRUIDA",
           area_ocupada = "AREA.OCUPADA",
           pavimentos = "QUANTIDADE.DE.PAVIMENTOS",
           ano_construcao = "ANO.DA.CONSTRUCAO.CORRIGIDO",
           tipo = "TIPO.DE.PADRAO.DA.CONSTRUCAO") %>% 
    
    # Separação do número de contribuinte (SQL) em setor quadra e lote
    mutate(setor =  str_sub(sql, 1, 3),
           quadra = str_sub(sql, 4, 6),
           
           # Quando o lote é um condomínio, haverá vários SQLs no mesmo lote. CD = Condomínio
           lote = str_sub(sql, 7, 10) %>% 
             ifelse(condominio == "00-0", ., paste("CD", str_sub(condominio, 1, 2), sep = "")),
           
           # Tipo de uso
           residencial = str_detect(tipo, "Residencial"),
           
           # Medida de verticalização
           verticalizacao = pavimentos * area_ocupada / area_terreno)
sql condominio area_terreno area_construida area_ocupada pavimentos ano_construcao tipo setor quadra lote residencial verticalizacao
0440880046-5 01-9 3760 211 1022 32 2013 Residencial vertical - padrão D 044 088 CD01 TRUE 8.6978723
1404360004-5 00-0 110 80 40 2 1997 Residencial horizontal - padrão B 140 436 0004 TRUE 0.7272727
0151260072-5 00-0 588 450 312 2 1989 Residencial horizontal - padrão D 015 126 0072 TRUE 1.0612245
1360660013-4 00-0 250 143 127 1 1991 Residencial horizontal - padrão B 136 066 0013 TRUE 0.5080000
0701830015-0 00-0 117 136 78 2 1977 Residencial horizontal - padrão C 070 183 0015 TRUE 1.3333333
0151020023-1 00-0 300 120 60 2 1957 Residencial horizontal - padrão C 015 102 0023 TRUE 0.4000000
0801210119-1 02-7 700 45 500 5 1970 Residencial vertical - padrão B 080 121 CD02 TRUE 3.5714286
1201400021-1 00-0 256 227 122 2 1978 Residencial horizontal - padrão C 120 140 0021 TRUE 0.9531250

Tipos de uso para IPTU

Todos os empreendimentos não públicos regulares (de acordo com a lei) constam nesta base de dados. A seguir está uma análise de quais são estes tipos. O tipo de interesse para esta análise é apenas os residenciais.

gg.right <- IPTU %>% 
  mutate(uso = case_when(str_detect(tipo, "Residencial") ~ "Residencial",
                         str_detect(tipo, "Comercial") ~ "Comercial",
                         str_detect(tipo, "Oficina") ~ "Industria",
                         str_detect(tipo, "TERRENO") ~ "Terreno",
                         str_detect(tipo, "Clube") ~ "Entretenimento",
                         .default = "Outros") %>% as.factor(),
         padrao = case_when(str_detect(tipo, "A$") ~ "A",
                            str_detect(tipo, "B$") ~ "B",
                            str_detect(tipo, "C$") ~ "C",
                            str_detect(tipo, "D$") ~ "D",
                            str_detect(tipo, "E$") ~ "E",
                            .default = "NA"),
         nome = paste("(", padrao, ") ", uso, sep = "")) %>% 
  group_by(uso, padrao, nome) %>% 
  summarize(area = sum(area_construida)) %>% 
  ungroup() %>% 
  mutate(percentual = area * 100/ sum(area, na.rm = TRUE),
         texto = case_when(percentual < 1 ~ "<1%",
                           percentual < 5 ~ paste(round(percentual,1), "%", sep = ""),
                           .default = paste(round(percentual, 2), "%", sep = "")),
         color = ifelse(uso == "Outros", "white", "black")) %>% 
  ggplot() +
  treemapify::geom_treemap(aes(fill = uso, area = area), size = 2, color = NA) +
  treemapify::geom_treemap(aes(alpha = padrao, area = area), fill = "black", color = NA) +
  treemapify::geom_treemap_text(aes(area = area, label = nome, color = color), alpha = .4, grow = FALSE, size = 8) +
  treemapify::geom_treemap_text(aes(area = area, label = texto, color = color), 
                                alpha = .4, place = "middle", size = 10) +
  scale_fill_manual(values = c("Residencial" = "#FAE48B", 
                               "Comercial" = "#A6EEE6",
                               "Industria" = "#84BD00",
                               "Entretenimento" = "#FB6467",
                               "Outros" = "#1A5354")) +
  scale_alpha_manual(values = c("A" = 0, "B" = .05, "C" = .1, "D" = .15, "E" = .2, "NA" = 0)) +
  scale_color_manual(values = c("black" = "black", "white" = "white")) +
  labs(fill = "Tipo de uso", alpha = "Parão de uso") +
  theme(aspect.ratio=1)

gg.left <- IPTU %>% 
  mutate(uso = case_when(str_detect(tipo, "Residencial") ~ "Residencial",
                         str_detect(tipo, "Comercial") ~ "Comercial",
                         str_detect(tipo, "Oficina") ~ "Industria",
                         str_detect(tipo, "TERRENO") ~ "Terreno",
                         str_detect(tipo, "Clube") ~ "Entretenimento",
                         .default = "Outros") %>% as.factor()) %>% 
  group_by(uso) %>% 
  summarize(area = sum(area_construida, na.rm = TRUE)) %>% 
  ungroup() %>% 
  arrange(desc(uso)) %>% 
  mutate(percentual = area * 100 / sum(area, na.rm = TRUE),
         texto = ifelse(percentual > 5, paste(round(percentual, 1), "%", sep = ""), ""),
         ytext = (cumsum(area) + lag(cumsum(area)))/ 2) %>% 
  ggplot() +
  geom_col(aes(y = area, fill = uso, x = "")) +
  geom_text(aes(y = ytext, fill = uso, x = "", label = texto), alpha = .5) +
  scale_fill_manual(values = c("Residencial" = "#FAE48B", 
                               "Comercial" = "#A6EEE6",
                               "Industria" = "#84BD00",
                               "Entretenimento" = "#FB6467",
                               "Outros" = "#1A5354")) +
  theme_void() +
  theme(legend.position = "none")

ggsave("tex/imagens/tree_area_construida.png", 
       (gg.left | gg.right) + 
         plot_layout(widths = c(1,6)), 
       width = 10, height = 8, dpi = 250)

Agrupamento dos condomínios

Quando há um condomínio com múltiplos contribuintes de IPTU, ele deve ser agregado a nível lote, para que possa ser cruzado com a base de lotes, possibilitando que seja georreferenciada. Todos os contribuintes dentro de um condomínio compartilham das mesmas características de área de terreno, área ocupada, ano de construção e pavimentos, então a mediana funciona para agregar estes dados, mas o máximo, mínimo, média ou pegar o primeiro valor também funcionaria.

IPTU <- IPTU %>% 
  group_by(setor, quadra, lote) %>% 
  
  # Agregar por SQL
  group_by(setor, quadra, lote) %>% 
  summarize(unidades = n(),
            area_terreno = median(area_terreno), 
            area_construida = sum(area_construida), 
            area_ocupada = median(area_ocupada),
            pavimentos = median(pavimentos),
            ano_construcao = median(ano_construcao),
            residencial = median(residencial),
            verticalizacao = median(verticalizacao)) %>% 
  ungroup()
setor quadra lote unidades area_terreno area_construida area_ocupada pavimentos ano_construcao residencial verticalizacao
029 022 0043 1 104 180 60 3 1954 1 1.7307692
055 179 0011 1 640 0 0 0 0 0 0.0000000
107 064 0008 1 200 225 113 2 1999 1 1.1300000
112 425 0049 1 378 384 378 1 1985 0 1.0000000
019 043 0036 1 90 100 50 2 1970 1 1.1111111
066 149 0088 1 96 47 47 1 1973 1 0.4895833
118 342 0041 1 480 250 150 2 2010 1 0.6250000
119 351 0023 1 135 195 96 2 2000 1 1.4222222
IPTU %>% 
  mutate(unitario = unidades == 1) %>% 
  group_by(unitario) %>% 
  summarize(n = n(),
            soma = sum(unidades))
## # A tibble: 2 × 3
##   unitario       n    soma
##   <lgl>      <int>   <int>
## 1 FALSE      33129 1782366
## 2 TRUE     1314353 1314353
gg <- IPTU %>% 
  ungroup() %>% 
  mutate(CA = area_construida / area_terreno,
         cota_parte = area_terreno / unidades) %>% 
  filter(cota_parte < 600, CA < 5, verticalizacao < 5) %>% 
  select("Cota Parte" = cota_parte, "Coeficiente de Aproveitamento (CA)" = CA, "Verticalização" = verticalizacao) %>% 
  pivot_longer(everything()) %>% 
  ggplot() +
  geom_histogram(aes(x = value)) +
  facet_wrap(~ name, scales = "free_x") +
  labs(title = "Distribuição dos indicadores em SP", x = "", y = "Frequência") + 
  theme_classic()

ggsave("tex/imagens/indicadores.png", gg, width = 9, height = 4, dpi = 250)

Geosampa

GeoSampa é o portal de mapas e dados geoespaciais da cidade de São Paulo, mantido pela prefeitura. Ele fornece uma vasta gama de informações geográficas, incluindo mapas, dados demográficos, infraestruturas e muito mais. Este portal é uma ferramenta valiosa para pesquisadores, urbanistas e qualquer pessoa interessada em informações espaciais detalhadas sobre a cidade.

https://geosampa.prefeitura.sp.gov.br/

Lotes

Base de Lotes: No GeoSampa, a base de lotes representa a divisão da cidade em pequenos segmentos, geralmente correspondentes a terrenos individuais ou condomínios. Esta base está organizada de forma que cada bairro tem seu próprio conjunto de dados e os dados de lotes para cada bairro podem ser baixados diretamente do site do GeoSampa em formato zip, contendo arquivos como .shp (shapefile), .dbf (database file), e .shx (index file).

# Extração do zip file com lotes de cada bairro
if (!"zip" %in% list.files(path = "dados/lotes")){
  for (file in list.files(path="dados/lotes/zip", full.names = FALSE) %>% 
       str_remove("\\.zip")){
    unzip(paste("dados/lotes/zip/", file, ".zip", sep = ""), 
          paste(file, "/", file, ".shp", sep = ""), 
          exdir = "dados/lotes/unzip")
    unzip(paste("dados/lotes/zip/", file, ".zip", sep = ""), 
          paste(file, "/", file, ".dbf", sep = ""), 
          exdir = "dados/lotes/unzip")
    unzip(paste("dados/lotes/zip/", file, ".zip", sep = ""), 
          paste(file, "/", file, ".shx", sep = ""), 
          exdir = "dados/lotes/unzip")
  }
}
# Junção dos shapes dos lotes de cada bairro em uma tabela
lotes <- list.files(path="dados/lotes/unzip", full.names = FALSE) %>% 
    paste("dados/lotes/unzip/", ., "/", ., ".shp", sep = "") %>% 
    lapply(read_sf) %>% 
    bind_rows %>% 
    st_set_crs("epsg:31983") 
lo_setor lo_quadra lo_lote lo_condomi lo_tp_quad lo_tp_lote geometry
149 057 0002 00 F F POLYGON ((347333.9 7390815,…
159 170 0004 00 F F POLYGON ((321624.4 7391258,…
085 001 0022 00 F F POLYGON ((327596 7386024, 3…
091 204 0045 00 F F POLYGON ((331395.1 7382865,…
195 007 0050 00 F F POLYGON ((352309.7 7385411,…
070 199 0003 00 F F POLYGON ((334931.9 7402747,…
149 305 0031 00 F F POLYGON ((346689.8 7390317,…
051 155 0044 00 F F POLYGON ((341296.9 7389320,…

Tipos de lotes

Os lotes podem ser classificados de três formas

  • Lotes fiscais: apresentam contribuintes do IPTU
  • Espaço livre: espaços públicos
  • Via de acesso: Ruas fechadas para acesso dentro de uma quadra

Os lotes fiscais podem ser unidades ou condomínios. Um prédio, por exemplo, fica em um único lote, mas dentro pode haver diversas unidades, então se configura como um condomínio. Não é possível saber através dos dados de lotes quantas unidades estão no condomínio, nem alguma forma de discriminá-los.

lotes %>% 
  mutate(condominio = case_when(lo_condomi == "00" ~ "Unidade",
                                lo_condomi != "00" ~ "Condomínio"),
         tipo = case_when(lo_tp_lote == "F" ~ "Fiscal",
                          lo_tp_lote == "M" ~ "Espaço livre",
                          lo_tp_lote == "V" ~ "Via de acesso",
                          .default = NA),
         area = st_area(geometry) %>% as.numeric()) %>% 
  st_drop_geometry() %>% 
  ggplot() +
  geom_violin(aes(x = factor(condominio), y = area, fill = tipo)) +
  scale_y_log10(labels = scales::comma_format(big.mark = ".")) +
  labs(title = "Área dos lotes em SP" , x = "", y = "Área em metros quadrados", fill = "Tipo de lote") +
  theme_classic()

SQL dos Lotes (Setor, Quadra e Lote)

Na base de dados de lotes, cada lote é identificado por três componentes principais:

  1. Setor (lo_setor): Uma divisão maior dentro do município que agrupa várias quadras.
  2. Quadra (lo_quadra): Uma subdivisão dentro de um setor que agrupa vários lotes.
  3. Lote (lo_lote): A menor unidade de divisão, que representa um terreno ou uma parcela específica dentro de uma quadra. Essa estrutura de Setor-Quadra-Lote é crucial para identificar de forma única cada lote dentro do município. No código, esses identificadores são utilizados para manipular e cruzar os dados de lotes com outras bases de dados, como a base de IPTU.
lotes <- lotes %>% 
 filter(lo_tp_lote == "F") %>% # Seleção apenas de lotes fiscais
  mutate(lo_lote = ifelse(lo_lote == "0000", paste("CD", lo_condomi, sep = ""), lo_lote)) %>% 
  select(setor = lo_setor, quadra = lo_quadra, lote = lo_lote)
setor quadra lote geometry
139 466 0015 POLYGON ((354220.7 7398122,…
175 236 0054 POLYGON ((327682.3 7374400,…
235 088 0036 POLYGON ((357544.1 7392499,…
016 162 CD04 POLYGON ((329068.7 7390030,…
114 328 0210 POLYGON ((352906.6 7397697,…
070 418 0032 POLYGON ((334452 7402651, 3…
089 465 0080 POLYGON ((332208.6 7385030,…
045 100 CD04 POLYGON ((332530.8 7387987,…

Base de quadras

quadras <- read_sf("dados/quadras/SIRGAS_SHP_quadraMDSF.shp") %>% 
  st_set_crs("epsg:31983") %>% 
  filter(qd_tipo == "F") %>% # Apenas quadras fiscais
  select(setor = qd_setor, quadra = qd_fiscal) %>% 
  group_by(setor, quadra) %>% 
  summarize(geometry = st_union(geometry)) %>% 
  ungroup()

Base de setores

setores <- read_sf("dados/setor/SIRGAS_SHP_setorfiscal.shp") %>% 
  st_set_crs("epsg:31983") %>% 
  select(setor = st_codigo) %>% 
  group_by(setor) %>% 
  summarize(geometry = st_union(geometry)) %>% 
  ungroup()

Join entre a Base de Lotes e a Base do IPTU

# Join dos lotes com IPTU com base no SQL
IPTU.lote <- IPTU %>% 
  filter(residencial == 1) %>% 
  left_join(lotes, by = join_by(setor, quadra, lote)) %>% 
  ungroup()
setor quadra lote unidades area_terreno area_construida area_ocupada pavimentos ano_construcao residencial verticalizacao geometry
152 286 0028 1 144 60 60 1 1984 1 0.4166667 POLYGON ((351034.9 7386395,…
115 310 0161 1 262 75 75 1 1970 1 0.2862595 POLYGON ((354872.7 7395381,…
085 021 0040 1 214 170 105 2 1973 1 0.9813084 POLYGON ((327642.7 7386172,…
152 391 0044 1 260 723 241 3 2015 1 2.7807692 POLYGON ((346753.2 7387844,…
166 155 0071 1 259 86 86 1 1980 1 0.3320463 POLYGON ((319217.7 7380822,…
057 045 0002 1 188 129 129 1 1976 1 0.6861702 POLYGON ((343003.4 7396046,…
065 078 0047 1 179 160 156 1 1954 1 0.8715084 POLYGON ((337017.7 7399615,…
049 048 0028 1 324 230 230 1 1970 1 0.7098765 POLYGON ((336888.9 7386628,…
IPTU.quadra <- IPTU %>% 
  filter(residencial == 1) %>% 
  group_by(setor, quadra) %>% 
  summarize(unidades = sum(unidades),
            area_construida = sum(area_construida),
            verticalizacao = weighted.mean(verticalizacao, area_terreno),
            area_terreno = sum(area_terreno)) %>% 
  left_join(quadras, by = join_by(setor, quadra)) %>% 
  ungroup()
IPTU.setor <- IPTU %>% 
  filter(residencial == 1) %>% 
  group_by(setor) %>% 
  summarize(unidades = sum(unidades),
            area_construida = sum(area_construida),
            verticalizacao = weighted.mean(verticalizacao, area_terreno),
            area_terreno = sum(area_terreno)) %>% 
  left_join(setores, by = join_by(setor)) %>% 
  ungroup()
list(IPTU %>% filter(residencial == 1) %>% mutate(base = "bruta"), 
     IPTU.setor  %>% mutate(base = "Setor")  %>% filter(!geometry %>% st_is_empty()), 
     IPTU.quadra %>% mutate(base = "Quadra") %>% filter(!geometry %>% st_is_empty()), 
     IPTU.lote   %>% mutate(base = "Lote")   %>% filter(!geometry %>% st_is_empty())) %>% 
  lapply(function(x) x %>% 
           group_by(base) %>% 
           summarize(unidades = sum(unidades))) %>% 
  bind_rows %>% 
  pivot_wider(names_from = base, values_from = unidades) %>% 
  pivot_longer(2:4) %>% 
  mutate(missing = bruta - value,
         missing_percent = (100 *missing / bruta) %>% round(2) %>% paste("%")) %>% 
  select("Critério de join" = name,
         "Erros (unidades)" = missing,
         "Percentual de erro" = missing_percent)
## # A tibble: 3 × 3
##   `Critério de join` `Erros (unidades)` `Percentual de erro`
##   <chr>                           <int> <chr>               
## 1 Setor                           -6409 -0.24 %             
## 2 Quadra                       -1025299 -38.56 %            
## 3 Lote                            44319 1.67 %

Algumas estatísticas descritivas

gg <- IPTU.lote %>% 
  mutate(distancia_centro = st_distance(geometry, 
                                        read_sf("dados/distrito/SIRGAS_SHP_distrito.shp") %>% 
                                          st_set_crs("epsg:31983") %>% 
                                          filter(ds_nome == "SE")) %>% 
           as.numeric() %>% 
           cut(breaks = 10^3*(0:40), labels = FALSE)) %>% 
  st_drop_geometry() %>% 
  group_by(distancia_centro) %>% 
  summarize(verticalizacao = weighted.mean(verticalizacao, area_terreno)) %>% 
  ggplot(aes(x = distancia_centro, y = verticalizacao)) +
  geom_col() +
  geom_hline(yintercept = 1, linetype = "dotted") +
  labs(x = "Distância ao centro em quilômetros", 
       y = "Verticalização") + 
  scale_y_continuous(breaks = (1:6)) +
  theme_classic()

ggsave("tex/imagens/verticalizacao.png", gg, width = 6, height = 4, dpi = 250)

gg.lotes <- IPTU.quadra %>% 
  ggplot() +
  geom_sf(data = distrito, fill = "black") +
  geom_sf(aes(fill = verticalizacao, geometry = geometry), color = NA) +
  theme_void() +
  scale_fill_gradient(low = "#004D40FF", high = "#E0F2F1FF")

gg.setor <- IPTU.setor %>% 
  ggplot() +
  geom_sf(data = distrito, fill = "black") +
  geom_sf(aes(fill = verticalizacao, geometry = geometry), color = NA) +
  theme_void() +
  scale_fill_gradient(low = "#004D40FF", high = "#E0F2F1FF")

ggsave("tex/imagens/mapa_verticalizacao.png", 
       (gg.lotes | gg.setor), width = 20, height = 20, dpi = 300)

Cruzamento entre a Base do Censo e a Base do IPTU

O processo de cruzamento foi realizado com base na intersecção das geometrias dos setores censitários e dos lotes do IPTU. Cada setor censitário e cada lote do IPTU possui uma geometria associada, representando sua área geográfica no mapa. Ao cruzar essas geometrias, é possível identificar quais lotes estão contidos em cada setor censitário e vice-versa.

É importante destacar que, em casos onde um lote foi dividido entre dois ou mais setores censitários, ocorrerá uma intersecção em ambas as áreas. Para lidar com essa situação, foi calculado o percentual da área do lote que está contida em cada setor censitário.

# Join dados do IPTU com do Censo através da intersecção das geometrias
IPTU.censo <- censo %>% 
  st_intersection(IPTU.lote %>% st_as_sf()) %>% 
  as_tibble() %>% 
  rename(geometria_intersec = geometry) %>% 
  
  # Retomada das geometrias do setor e lotes
  left_join(censo %>% 
              as_tibble() %>% 
              select(id_setor_censitario, 
                     geometria_setor_censitario = geometry),
            by = join_by(id_setor_censitario)) %>% 
  left_join(IPTU.lote %>% 
              as_tibble() %>% 
              select(setor, quadra, lote, 
                     geometria_lote = geometry),
            by = join_by(setor, quadra, lote)) %>% 
  
  # Cálculo de quanto % do lote está dentro do setor
  mutate(percent_intersec = as.numeric(st_area(geometria_intersec) / st_area(geometria_lote))) 
id_setor_censitario v0001 v0002 v0003 v0004 v0005 v0006 v0007 area_setor setor quadra lote unidades area_terreno area_construida area_ocupada pavimentos ano_construcao residencial verticalizacao geometria_intersec geometria_setor_censitario geometria_lote percent_intersec
355030893000182P 650 313 313 0 2.407407 3.703704 270 66545.90 [m^2] 044 092 0018 1 360 600 300 2 2019 1 1.6666667 POLYGON ((338893.1 7390604,… MULTIPOLYGON (((338978 7390… POLYGON ((338891.5 7390590,… 1.0000000
355030896000048P 765 271 270 1 3.031873 5.179283 251 47687.19 [m^2] 115 343 0024 1 115 200 100 2 1991 1 1.7391304 POLYGON ((356750.5 7396784,… MULTIPOLYGON (((356861.6 73… POLYGON ((356749.5 7396782,… 1.0000000
355030828000066P 728 323 323 0 2.545455 8.391608 286 34056.15 [m^2] 130 286 0014 1 299 400 210 2 1991 1 1.4046823 POLYGON ((347710 7401170, 3… MULTIPOLYGON (((347803.1 74… POLYGON ((347713.2 7401171,… 1.0000000
355030871000023P 384 164 164 0 2.844444 28.888889 135 106275.40 [m^2] 088 265 0351 1 401 410 205 1 1986 1 0.5112219 POLYGON ((327758.4 7384904,… MULTIPOLYGON (((327900 7385… POLYGON ((327754 7384919, 3… 1.0000000
355030896000071P 608 234 234 0 2.881517 6.161137 211 26319.25 [m^2] 138 176 0026 1 125 119 88 2 2001 1 1.4080000 POLYGON ((355278.8 7397050,… MULTIPOLYGON (((355288.8 73… POLYGON ((355277.1 7397045,… 1.0000000
355030868000149P 748 303 303 0 2.888031 9.266409 259 25093.66 [m^2] 157 020 0035 1 140 143 79 2 1982 1 1.1285714 POLYGON ((335749.3 7384182,… MULTIPOLYGON (((335769.3 73… POLYGON ((335757.4 7384187,… 1.0000000
355030864000080P 548 212 212 0 2.795918 7.142857 196 30000.23 [m^2] 111 151 0002 1 220 75 50 2 1985 1 0.4545455 POLYGON ((348713.8 7398937,… MULTIPOLYGON (((348791.1 73… POLYGON ((348735.6 7398930,… 1.0000000
355030886000085P 581 271 271 0 2.702326 6.976744 215 82078.28 [m^2] 068 555 0017 1 301 90 90 1 1970 1 0.2990033 POLYGON ((335447 7400688, 3… MULTIPOLYGON (((335682.5 74… POLYGON ((335448.7 7400679,… 1.0000000
355030821000028P 586 312 312 0 2.700461 8.294931 217 54365.98 [m^2] 075 258 0044 1 252 332 95 3 2002 1 1.1309524 POLYGON ((331400.4 7400874,… MULTIPOLYGON (((331685.8 74… POLYGON ((331398.2 7400876,… 1.0000000
355030829000026P 266 115 115 0 2.891304 32.608696 92 21554.58 [m^2] 104 077 0026 1 300 80 80 1 1988 1 0.2666667 POLYGON ((326127.3 7400794,… MULTIPOLYGON (((326131.9 74… POLYGON ((326129.8 7400790,… 0.9884079

Erros no join

Disparidades de informações

Teoricamente o número de domicílios deveria ser próximo do número de unidades habitacionais enunciadas pelo IPTU. Caso haja mais unidades habitacionais do que domicílios, houve algum equívoco na medição do censo, visto que entre domicílios não contam apenas os ocupados. Caso haja mais domicílios, há unidades que não são contribuintes do IPTU. O que é preocupante é que casos em que o número é diferente não são exceções.

gg <- censo %>% 
  st_drop_geometry() %>% 
  select(id_setor_censitario, domicilios = v0002) %>% 
  left_join(IPTU.censo %>% 
              st_drop_geometry() %>% 
              filter(residencial == 1) %>% 
              group_by(id_setor_censitario) %>% 
              summarize(unidades = sum(unidades))) %>% 
  mutate(unidades = replace_na(unidades, 0)) %>% 
  mutate(balanco = unidades / (unidades + domicilios)) %>% 
  ggplot() +
  geom_histogram(aes(x = balanco)) +
  labs(caption = "Caso o balanço seja 50%, implica que o número de domicílios reportados no censo é igual ao número de 
unidades reportadas no IPTU. Caso o balanço seja menor que 50%, implica que há mais domicílios no censo 
do que unidades, e caso seja maior que 50%, há mais unidades que domicílios",
title = "Disparidade entre censo e IPTU quanto ao número de unidades habitacionais",
    
       x = "Balanço: unidades / (unidades + domicilios)", y = "") +
  theme_classic() +
    theme(plot.caption = element_text(hjust = 0)) +
  scale_x_continuous(labels = scales::percent)


ggsave("tex/imagens/disparidade_censoIPTU.png", gg, width = 8, height = 5, dpi = 400)

Perfil geográfico dos erros

Alguns setores censitários não encontram pares na base de lotes. Isso acontece principalmente por conta de loteamentos irregulares e favelas, que não são contribuintes do IPTU, e, portanto, não constam na base. Isso não prejudica a análise, pois estes loteamentos não dependem da regulamentação urbana, então mudanças nos instrumentos e indicadores não impactariam essas regiões.

Outras falhas decorrem do erro do join da base do IPTU com lotes, como apontado anteriormente, mas estes casos são negligenciáveis.

distrito <- read_sf("dados/distrito/SIRGAS_SHP_distrito.shp") %>% st_set_crs("epsg:31983")
lotes_irregulares <- read_sf("dados/lotes_irregulares/SIRGAS_SHP_loteamento.shp") %>%
  st_set_crs("epsg:31983")
favela <- read_sf("dados/favela/SIRGAS_SHP_favela.shp") %>% 
  st_set_crs("epsg:31983")

censo.pontos <- censo %>%
  select(id_setor_censitario, geometry, pessoas = v0001, domicilios = v0002) %>%
  mutate(pontos = round(pessoas / 100)) %>% 
  as_tibble() %>% 
  left_join(
            # Seleção dos setores censitários que não apresentam nenhum contribuinte do IPTU
            censo %>% 
              anti_join(IPTU.censo) %>% 
              st_drop_geometry() %>% 
              select(id_setor_censitario) %>% 
              mutate(erro = TRUE)) %>% 
  mutate(erro = replace_na(erro, FALSE)) %>% 
  mutate(samples = map2(geometry, pontos, ~st_sample(.x, size = .y))) %>%
  unnest(cols = samples) %>%
  as_tibble()

gg <- censo.pontos %>% 
  ggplot() +
  geom_sf(data = distrito, color = NA) +
  geom_sf(data = st_union(favela %>% st_union() %>% st_buffer(10),
                          lotes_irregulares %>% st_union() %>% st_buffer(10)),
          aes(fill = "Favelas e lotes irregulares"), 
        color = "black", size = .1, alpha = .7) +
  geom_sf(aes(geometry = samples, color = erro), alpha = .25, size = .2) +
  scale_color_manual(values = c("FALSE" = NA,#248232
                               "TRUE" = "#D32934FF"),
                     labels = NULL) +
  scale_fill_manual("", values = c("Favelas e lotes irregulares" = "#2BAA92FF")) +
  labs(title = "<span style='font-size: 35pt;'>População em <span style = 'color:#D32934FF;'>áreas sem registro de IPTU</span> geralmente estão em <br><span style = 'color:#2BAA92FF;'>favelas ou lotes irregulares</span> (cada ponto representa 100 pessoas)</span>") + 
  theme_void() +
  theme(plot.title = element_markdown(), legend.position = "none")

# scale_colour_paletteer_d("lisa::AndyWarhol_2")

ggsave("tex/imagens/mapa_pontos.png", gg, width = 18, height = 20, dpi = 300)